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O , Abstract 

• ^ ' The steady-state simplified P/v approximation to the radiative transport equation 

has been successfully applied to many problems involving radiation. Recently, time- 
I dependent simplified Pjy equations have been derived by an asymptotic analysis 



Ph. 



similar to the asymptotic derivation of the steady-state SPn equations [7J. In this 
paper, we present computational results for the time-dependent SPn equations in 
two dimensions, obtained by using an adaptive finite element approach. Several 
numerical comparisons with other existing models are shown. 

in 
p 

O ; 1 Introduction 
a^ 
o 

Time-dependent radiative transfer, described by the radiative transfer equation, is hard to 
^ , compute. This is due to the six-dimensional phase space (Ix time, 2x angle, 3x space). 
^ I There is an interest in time-dependent radiative transfer solutions, e.g. in astrophysics 
(supernova explosions), the interaction of short-pulsed lasers with plasmas, and light de- 
tection and ranging (LIDAR). It is our purpose in this paper to numerically investigate 
the time-dependent simplified Pn equations in two dimensions. These models have been 
very successful in the steady case. Here we investigate an extension to the time- dependent 
case. 
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The simplified P/v (SPn) equations were originally developed for steady-state problems 
in nuclear engineering [H [9l [10] and have subsequently been generahzed and successfully 
applied in several other fields, including radiative transfer [13l [TBI [19]. The first formal 
derivation by Gelbard P, IHl HD] started with the one- dimensional P/v equations, which 
contain only first-order space derivatives, and used substitutions to obtain a system of 
elliptic partial differential equations. To obtain equations in three space dimensions, even- 
order moments are interpreted as scalars, odd-order moments are interpreted as vectors, 
and one-dimensional derivatives dx are replaced by divergence operators and gradients 
respectively. In three space dimensions, compared to the (A^ + 1)^ independent unknowns 
in the spherical harmonics P/v equations, the number of unknowns in the SPn equations 
increases only linearly with A^. Because of the derivation via the one-dimensional P/v 
equations, the SPjy method was at first not widely accepted. But alternative derivations 
via asymptotic expansion [17] and via a variational approach [21 [23] have substantiated the 
validity of the SP^ hierarchy. 

The SPn equations are accurate if the medium is optically thick, the scattering rate is 
comparable to the collision rate, and scattering is not highly forward-peaked [17]. In 
addition, numerical experiments (cf. ^19] and references therein) have shown that the SPn 
equations give good results even when the regime is not so diffusive, and even in the 
presence of a discontinuity in the opacities. This means that in the diffusive regime a 
higher accuracy is obtained and at the same time the range of applicability is increased. 

Until recently, the SPn method was almost exclusively applied to steady-state transport 
equations, i.e. no time dependence was assumed. Only then can the P/v equations be 
substituted into each other to give a second-order system. To our knowledge, there has 
been only one attempt in the literature [21j to apply the SPn method to a time- dependent 
problem. Here, the authors use a semi-discretization in time (i.e. the time variable is 
discretized whereas the other variables are treated as being continuous) and apply the SPn 
approximation to the then steady system. This paper, on the other hand, investigates time- 
dependent SPn equations which were systematically derived from the Boltzmann equation 
using an asymptotic analysis. 

The numerical solution of the SPn equations are often still quite expensive due to their 
inherent multi-scale structure in both time and space. A remedy is to use fully adaptive 
algorithms where the local accuracy of the numerical solution is controlled by means of a 
posteriori error estimates in space and time. Such estimators are well established to control 
the adaptive multilevel process producing highly refined space-time grids to capture local 
effects efficiently and therefore drastically reducing the size of the arising linear algebraic 
systems with respect to a prescribed tolerance. We apply the adaptive Rothe method 
based on the discretization sequence first in time then in space, in contrast to the usual 
Method of Lines approach (see e.g. [H] and references therein). The spatial discretization 
is considered as a perturbation of the time integration process. Implementations have been 
done in the KARDOS library [6], which provides a suitable programming environment for 
adaptive algorithms to solve nonlinear time-dependent PDEs. 
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This paper is organized as follows: A brief summary of the derivation of the time- dependent 
SPn equations using asymptotic analysis is given in Section [2l Suitable initial and bound- 
ary conditions are stated in Section [3l The imployed numerical method is described in 
Section HI In Section [5l these techniques are applied to two test cases from the recent 
radiative transfer literature. 



2 Time-Dependent SPn Equations 

We consider a convex, open, bounded domain Z in M^, and we assume that Z has a 
smooth boundary with outward normal vector n. The direction of particle motion is given 
by G 5*^, where 5*^ is the unit sphere in three dimensions. Moreover, we let 

T = dZ X and = {{x, ^]) G T : n{x) ■ < 0}. 

The transport of mono-energetic particles that undergo isotropic scattering in a medium 
is modeled by the linear Boltzmann equation 

-dttpit, x,n) + n- Va^tpit, X, n) + at{x)ij{t, x, n) = [ ^(t, x, n')dn' + ^i^, (2.1) 

V An J g2 An 

where q is an isotropic source term. At the boundary, we prescribe the ingoing radiation 

^{t,x,Vt) =^i,{t,x,VL) on r-, (2.2) 

and as the initial condition, we prescribe 

V'(0,a;,n) = ^/'o(x,fi). (2.3) 

Here, ip{t,x,VL) cos 9 dAdtdVL is the number of particles at point x and time t that move 
with velocity v during dt through an area dA into a solid angle dVL around f2, and 9 is the 
angle between VL and dA. The total cross section at{x) is the sum of the absorption cross 
section o"a(x) and the total scattering cross section (Js^x). 

The time-dependent SP^ equations have been derived in |7j. For the convenience of the 
reader, we here present an abbreviatied version which contains the major ideas. The steady- 
state diffusion equation is an elliptic PDE. Time- dependent diffusion theory is governed 
by a parabolic PDE. To obtain higher-order corrections to diffusion theory, we write the 
transport equation in a parabolic scaling. Space-derivatives are scaled by a small parameter 
e and the additional time-derivative is scaled by e^. This is called a parabolic scaling, since 
a differential operator that is first-order in time and second-order in space is invariant 
under this scaling. The transport equation is therefore written as: 

e^-dtij + en ■ + (Ttij = {at - eV,) -^0 + e^^, (2.4) 
V ^ ' An An 
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where ■?/' = ip(t, 2;, fl), (j)(t, x) = J^a "^{t, x, Q)dQ, and q = q(t, x). 
Integrating (12. 4p over fl and dividing by e"^, we obtain the "balance" equation 

1 1 f 

-dt(j) + -V:r- / QipdQ + aa4> = q, (2.5) 
V e Js2 

which states a basic physical principle: changes in the scalar flux are either due to 
leakage (the spatial derivative term), absorption, or sources. We require that this "balance" 
equation be contained in the final choice of SP^ equations. 

We write ([22]) as 

{l+en-X + £^T)iP = S, (2.6) 



where 



X = — V^, T= dt, and S = { 1 - e'— ] ^ + e'-^. (2.7) 

at vat V CTt J An Anat 



We start by expanding the inverse of the operator in fl2.6p in powers of e 
^ = {l + eQ-X + e^T)-^S 

= {i-{n-x)e+[-T + {n- xf] + [{n ■ x)t + {T-{n- xf){n ■ x)] 

+ [{T - {Q ■ Xf)T + (-2(fi • X)T + {Vl ■ Xf){VL ■X)\e^---]S + O(e^). (2.8) 

In the following we assume that the system is homogeneous, i.e. cTq and at are constant. 
This assumption is crucial for the validity of the following analysis. For a discussion of 
the non-homogeneous case, we refer the reader to the end of this section. Integrating (12. 8p 
with respect to VL and using 

/ {n . xrdn = [1 + (-1)"] X" = [1 + (-1)"] (X ■ x)t, (2.9) 

J 32 n + l n + 1 



we obtain 



ipdn 

S2 

=47r |l + Qx^ -T^e'+ (^T^ + ^X^ - TX^ ] 
+ (^X^ + 2T'^X^ -T'^- TXM e'^\S + 0{e^). (2.10) 
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Hence, 



-X^ - T]e^+\T^ + -X^ - TX' e 



+ ( -X*' + 2X^X2 - - TX" 1 e 



-1 



1 



--X^ + T\e^ + 



45 



X4 + ITX^ e 



+ 



44 



X^ 



4 



T'X' + — TX^ + C(e'). 



945 3 15 
Inserting the definition of tlie source term S from (12. 7p . we get 



.2^ 



1 



--X2 + T + -— + -TX^ e 



45 



+ 



44 



X« - -T^X^ + — TX^ 1 e 



945 3 15 
Deleting on both sides and muhiplying by cTt/e'^, we obtain 



-aa(j) + q = crtT(j) - yX^ 



15 



44 



315 



4 
5^ 



-e^TX^c 



(2.11) 



f2.12) 



(2.13) 



We note that this equation has the form of the balance equation (12.51) . Since we want to 
keep this form, in the subsequent approximations we only manipulate the terms within the 
brackets. 



2.1 SPi Approximation 

For the lowest-order approximation, we neg lect terms of order 0{e^). Then (12131) becomes 
the classical diffusion (SPi) equation 

-dt(l) = ^Vl<P-aa<P + q. (2.14) 



2.2 SPs Approximation 

As in the steady case, the SP2 equations, which are of order 0{e'^), have proven to be 
inadequate in practice. This is due to their origin from the P2 equations [Tj. Therefore we 
omit them and proceed with the SP3 equations. 
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Noting that Eq. fl2.13p has the form of the balance equation (12.51) . we write (I2.13P as 



q - CFaC, 



+ 



3ae^T 



15 



:e'X\ 



2 v2 



l-e'T + -(l-a)e'X 



(2.15) 



As before, we have isolated terms that contain time- dependent diffusion operators (first- 
order time and second-order space derivative). The transformation of the asymptotic 
expansion into the SP2 system, i.e. the definition of ^, is unique up to a multiplicative 
factor. However, for the expansion up to terms of order 0{e^), it is not clear how the 
substitutions have to be performed. Thus we have introduced a parameter a G [0, 1] to 
split the mixed term TX^ into two parts. We chose the parameter between zero and one 
in order to get diffusion equations with the correct signs. 



Using Neumann's series, we write (12.151) as: 



q-aa(j) = (ytT<f) - yX^ 



+ 



1 _ 11^2x2 + Sae'T 
21 

l + e^T--il-a)e''X^ 
5 



-1 4 

15" 

-1 



e'X' 



+0(56). 



Now we define 



1 - —e^X^ + 3ae^T 
21 



l + e'T--{l~a)e^X^ 
5 



15 ^ 



to obtain the system 



1 

3(Xt 

3at ^ 
V 3at 



Vl [(f) + 202 - C] - cra(j) + g. 



W2 



I5a 

+ ^ 



11 ^ 

21a 



3a £2 



12 

^2+ I — (1 -«) 



0-* 



(2.16) 

(2.17a) 
(2.17b) 

(2.18a) 
(2.18b) 
(2.18c) 



Without time- dependence, the variable C. is zero. Moreover, for a = | the above equations 
reduce to the steady-state SP3 approximation. To obtain a system that is not ill-posed, 
we must take < a < 0.9 [7|. 



2.3 Simplification of the SP^ System 

In [7j, the SP3 equations with a = 2/3 were derived from the P3 moment equations. 
The variable 02 can be identified with the second-order Legendre moment of the radiative 



6 



intensity. The variable (, on the other hand, is an auxihary variable without a straight- 
forward physical interpretation. Furthermore, C = in steady-state. To simplify the SP^ 
equations, we therefore make a quasi-steady approximation and neglect (. We obtain 

-9i0=-^V^[0 + 202]-a„0 + g, (2.19a) 
V Sat 



-dth = ^^l 
V Sat 



2 11 

+ 



3^>- (2.Wb) 



15a 21a 

We call these the SSP3 (simplified-simplified P3) equations. 

We expect that the time-dependent SP^ equations can be generalized to anisotropic scat- 
tering in a similar manner as in the steady-state case [17]. In the derivation of the equations, 
we assumed a homogeneous medium. In steady-state, a variational analysis yielded the 
SPn equations for non-homogeneous media as well as interface and boundary conditions 
[2], [23] . The only difference for space-dependent coefficients is that the spatial derivatives 
have to be modified like ^ ^ 

-v^ ^ v.— -v.. 

CTt crt{x) 

For steady-state problems, this modification of the spatial derivatives is asymptotically cor- 
rect in planar geometry and we expect that it is asymptotically correct for time-dependent 
planar geometry problems. 



3 Boundary Conditions and Initial Values 

In this section, we state boundary conditions for the SPn equations which have been 
derived using Marshak's method |20]. Let 



JnfKO JnQ<0 

For the SPi equations, we have 



e \2 2 

For SP3 and SSP3, we obtain the boundary conditions 



n-V.0=^fii-^0). (3.1) 



n ■ = ^ + + Ik + ^2) (3.2a) 

n-V.0.^^(l0-|0.-l/.) (3.2b) 
C = 0. (3.2c) 
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3.1 Initial Values 



Given an initial particle distribution, it is straight-forward to calculate an initial value 
for (f). From the asymptotic analysis, the physical meaning of the auxiliary variables 
02, C) is not obvious. Therefore it is not clear what the appropriate initial conditions 
for these variables are. In many cases, the initial setting is a steady state. In addition, 
the time-dependent SPjq equations reduce to the steady-state SPj^ equations. For the 
SP^ equations, we would have to solve (12.171) for 02 and C- Of course, this gives C = 0- 
Alternatively, 02 could be identified as the second-order Legendre moment and thus be 
computed from the initial value for ip. 

4 Numerical Method 

We have derived the time- dependent S'Pat equations in three spatial dimensions. In the 
following chapter, we will present numerical results in two spatial dimensions. Therefore, 
here we also describe the numerical method for two dimensions. It can be generalized, 
however, in a straightforward way. 

Mathematically speaking, the above SP^ models are nonlinear parabolic PDEs. This 
means that after spatial discretization we are faced with a large scale stiff system. From 
an ODE point of view, an optimal treatment of this stiffness structure is to apply some 
L-stable implicit time discretization [H [12]. From the PDE point of view, the avoidance 
of order reduction (which may occur above order 2) is equally important for the overall 
efficiency of the time integrator. Both properties are satisfied by the linearly implicit time 
discretization of Rosenbrock type behind the code ROSS PL [15]. Note that the popular 
Crank-Nicolson scheme is not L-stable and even not strongly A-stable, which results in an 
insufficient filtering of spurious modes. Fully implicit schemes require the iterative solution 
of finite dimensional nonlinear systems of algebraic equations by some Newton-like method. 
In contrast, linearly implicit methods realize a simplified Newton method in function space 
and require only the solution of a fixed number of linear systems per time step. According to 
their one-step nature, they allow for a rapid change of step sizes and an efficient adaptation 
of the spatial discretization in each time step. Moreover, a simple embedding technique 
can be used to estimate the error part arising from time discretization. 

Linearly implicit time integrators of Rosenbrock type are implemented in the code fam- 
ily KARDOS [6], which is used to solve our SP^ models. KARDOS is characterized by 
a combination of Rosenbrock solvers in time with multilevel finite elements in space in 
the setting of an adaptive Rothe approach, i.e., first time discretization and then spatial 
discretization. In this setting, both time-step control and dynamic mesh refinement on 
the basis of a posteriori error estimation can be simultaneously realized [H US]. A rigorous 
analysis for nonlinear parabolic systems has been given in [H], where challenging examples 
from other fields of science and technology are also included. 
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Next we want to describe the main ingredients of the adaptive Rothe method as needed 
for the efficient solution of the above described SP^ models. These models can be written 
as abstract Cauchy problems of the form 

HdtU = F{U), U{to)=Uo, t>0, (4.1) 

where if is a constant regular matrix and the diffusion operators and the boundary 
conditions are incorporated into the nonlinear function F{U). For example, we have 
U = (0, ^, 0'^ for the SP3 model. To approximate the vector U{x, t) by values Un ~ U{-,tn) 
at a certain time grid 

= to < ti < ■ ■ ■ < tn < ■ ■ ■ < tM~i <tM = T, (4.2) 

we apply the 4-stage third-order Rosenbrock method R0S3PL, which has the recursive 
form 

Jn]Um = F[Un + Y,('^lUnJ] -Hj2^Unj, 2 = 1, . . . , 4, (4.3) 





n 



Un+1 = t/„ + ^mi?7™, (4.4) 

i=l 

where r„ = — 1„ and J„ = F'{Un). The defining formula coefficients mj, aij, Cij, and 
7 are given in [15]. The method is L-stable and avoids order reduction. 

R0S3PL offers a simple way to estimate the local error. An embedded solution Un+i of 
second order can be computed by replacing the original weights rrii by rhi in (14.41) . In order 
to take into account the scale of the problem, the local error estimator is defined by the 
weighted root mean square norm 

( \\Un+l - Un+l\\h(^Z) \ ' , . 

\atol + rtol\\u,,^41,^,^) ■ ^^-^^ 

The tolerances ATOL and RTOL have to be selected carefully to furnish meaningful input 
for the error control. The estimator can be used to propose a new time step by 

Tn fTOLtrnV'^ , 

Tn+1 = Tn, (4.6) 

where TOLt is a desired tolerance prescribed by the user [TT]. If r„+i > TOLt, the step is 
rejected and redone. Otherwise the step is accepted and we advance in time. 

Observe that the above time discretization scheme has been applied to the abstract Cauchy 
problem (14.11) . i.e., to the initial value problem in function space. This means that the 
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Rosenbrock discretization scheme (14 .Sp is a sequence of linear elliptic boundary value prob- 
lems. The spatial approximation of the vectors f/„ is now done by multilevel finite elements. 
This is described next. 

The main idea of multilevel techniques consists of replacing the infinite dimensional solu- 
tion space by a nested sequence of finite dimensional spaces with successively increasing 
dimension in order to improve the approximation quality. To construct adaptive spatial 
meshes, we apply the edge-oriented hierarchical error estimator from [5l|ll]. Such estima- 
tors are well established to control the adaptive multilevel process producing successively 
finer meshes and, in spatial multi-scale cases, drastically reducing the size of the arising 
linear algebraic systems with respect to a prescribed tolerance. Let be an admissible 
finite element mesh at t = tn and S\ be the associated finite dimensional space consisting 
of all continuous piecewise linear functions. Then the standard Galerkin finite element 
approximation U!^^ G 5*^ of the intermediate values Uni in (14. 3 p satisfies the equation 

{Ln Ut. <Ph) = {Rni, (ph) for all 0,, G Sl (4.7) 

where L„ is the weak representation of the differential operator on the left-hand side in (14. 3 p 
and Rni stands for the entire right-hand side in (14. 3p . Since the operator L„ is independent 
of i, its calculation is required only once within each time step. The resulting large scale 
linear systems are solved by the BICGSTAB algorithm [24j with ILU preconditioning. 

After computing the approximate intermediate values f/^j, a posteriori error estimates 
can be used to give specific assessment of the error distribution. Consider a hierarchical 
decomposition 

Sl = Sl® Zl (4.8) 

where Z^ is the subspace that corresponds to the span of quadratic bubble functions 
corresponding to edges. Defining an a posteriori error estimator E^j^^ G by 

4 

i=l 

with approximating the projection error of the initial value Un in Z| and E^^ estimating 
the spatial error of the intermediate value ?7^j, the local spatial error for a finite element 
T G 7/i can be estimated by rjT '■= H-E^+iHt- The error estimator E'^j^^ is computed by 
linear systems which can be derived from (14.70 . We get for i = 

{LnEil„(j)h) = {Ln{Un-U!l),(j)H) forall0, gZ2. (4.10) 

and for i = 1, . . . , 4 

(L„ 0,) = iRr.iEl^^ + UH,, . . . , + t/i,_i), 0.) - (L. UH,, 0,) for all 0, G Zl 

(4.11) 
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Solving these equations encounters a sequence of five large linear problems in the space of 
hierarchical surpluses. From many practical computations, we have experienced that using 
the approximate error estimator 

E„V^Eti = E„^ + ^, (4.12) 

that is an error estimator for the embedded, locally second order linearly implicit Euler 
solution u^^^^^^"^ = -\- is quite efficient. Combining fl4.10p and the first equation 

of (14.111) yields the following simplified error equation 

(L„ E^+i, 0,) = (L„(t/„ - f/„\f ^0 + 0,) for all 0, G Zl (4.13) 

7 

Although we have reduced the number of error equations considerably, we still face a fully 
coupled system over the surplus space in (I4.13P . Following the approach given in [5] , we 
take further advantage of a localization strategy. The idea is to replace the bilinear form 
on the left hand side in f l4.13p by a spectrally equivalent block-diagonal preconditioner in 
the surplus space. Then, the error equation can be simultaneously solved for each bubble 
function, that is, for each edge in the triangulation. Let E^^^ be the corresponding error 
estimator. The local spatial error rjT for a finite element T &Th can again be estimated by 
computing the norm of E^'^^ over T. For the overall spatial error, we define in line with 
the local temporal error in (14.50 

\ATOL + RTOL\\UU\l^^z)) ' ^' ^ 

Based on this error estimation, we can control the spatial accuracy of the numerically 
computed solution to an imposed tolerance level TOL^. New grid points are placed in 
regions of insufficient accuracy. Therefore, all elements with rjT > CSmaxyr/j' are refined. 
We apply the standard red-green refinement technique. The iterative process estimate- 
refine-solve within a time step is continued until |||^^!^°^||| < TOL^. Obviously, temporal 
and spatial errors have to be well balanced. We have also to take into account mesh 
coarsening to gain efficiency. For more details, we refer to 



5 Numerical Results 
5.1 Marshak Wave 

This test case is a two-dimensional version of the analytical Marshak Wave test case from 
|22] . Here, radiation is coupled to an energy equation for i? ~ T^. The heat capacity is 
chosen such that the problem becomes linear. The equations are 
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(a) Cut along y = 0. (b) Cut along x = y. 

Figure 1: Energy distribution for different times in 2D Marshak wave. 



-dtij{t, x,Q)+Q- V^ij{t, X, Q) = aa{x) {B{t, x) - ij{t, x, Q)) + Q{t, x) (5.1) 
c 

dtB{t, x) = aa{x) (^-^ ipit, X, n)dn - B{t, x)^ . (5.2) 

The SPn approximation is applied to (15. ip and treats the B variable as an additional 
source term. The additional equation (15.21) is an ordinary differential equation but fits into 
the numerical framework above. 

The setting is two-dimensional and infinite in space {x G M^), with time t G [0, 10]. We 
have (Ta = 1.0. In an initially empty medium, a spatially bounded source Q is switched on 
at time zero: 



Q{t,x) = 




for < t < t*, X G [— xo,xo] x [— xo,a;o], 
otherwise 



with xo = 0.5 and t* = 10. 

For symmetry reasons, the computational domain was restricted to [0, 10] x [0, 10]. We set 
TOLt = TOL^ = 10-^ RTOL = 1, ATOL = 10"! The first time step was tq = l.Oe - 4. 
We started with a criss-cross grid of 20^ -|- 21^ = 841 points, the largest side length in one 
triangle being h = 2~^. After adaptive refinement, the smallest side length was h = 2~^°\/2- 
The finest grid consisted of 58916 points (both values extremal for SP3). Both spatial 
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(a) Mesh. (b) Time step. 

Figure 2: Spatial degrees of freedom (left) and step size (right) as functions of time. 



degrees of freedom and the time step are shown in Figure [21 The total computation times 
were 19m33s (SPi), 79ml0s (SSPg), 208m53s {SP3) on a PC with a 3 GHz 1686 processor. 

In Figure [H the radiative energy for times t = 1 (lower curves) and t = 10 (higher 
curves) computed by the different models 5* Pi, SP3, SSP^ is compared to a benchmark 
solution (high-order spherical harmonics solution), for both a cut along the well 
as a cut along the diagonal x = y. As in the one-dimensional case (investigated in [7], 
the higher-order diffusion approximations are a clear improvement on the SPi diffusion 
approximation. 



5.2 Lattice Problem 

As a more complex numerical test we consider a 2D checkerboard structure of different 
materials. The geometry, shown in Fig. [3], is identical to the example presented in [3], 
however modified here to have at = (Xg = 0.2 cm~^ in the highly scattering regions (white 
in Fig. [3]), and at = 10 cm~^, ag = cm~^ in the highly absorbing regions (grey and 



hatched squares in Fig. 4(a)). A source (hatched square in Fig. [3]) g = 1 is switched on 
at t = 0. The final time is t = 2.0. The test lies in an intermediate regime between thin 
and thick media. The radiation field propagates through several obstacles. The main front 
propagates to the top which is open, but some radiation will leak through the squares. 

We set TOLt = lO'^TOL^ = 10-^ RTOL = 1, ATOL = 10"^ with initial time step 
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TIME STEP SIZES 




(a) Mesh. (b) Time step. 

Figure 3: Computational mesh for the SP3 solution at t = 2.0 (left) and step size as a 
function of time (right). 
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Icm 1cm 1cm 1cm 1cm 1cm 1cm 

(a) Setting. (b) Cut along x = 3.5. 

Figure 4: Computational domain for the SP3 solution (left) and comparison of the different 
radiative energies along the center line x = 3.5. 
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(c) SSP3 solution. (d) SP3 solution. 

Figure 5: Energy distribution at t = 2.0 for different models. 
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To = 10~^. Figure 3(a) shows the adaptive grid for the SP3 solution at t = 2.0. It consists 
of 100100 triangles and 50135 points, the smallest triangle being 2~^. Local refinements 
can be seen around the jumps at the interfaces between the two media and in the central 
radiating region. Computation times were 23m46s (SPi), 58ml6s {SPP3), 164m57s {SP3) 
on a PC with a 3 GHz 1686 processor. 

Figure [5] shows the radiative energy at time t = 2.0 on a logarithmic scale. The P7 solution, 
which serves as a benchmark here, shows sharp fronts filling the squares adjacent to the 
center square and a small front escaping on the top. While the radiation in the center region 
is sufficiently well-described by the diffusion solutions, they overpredict the spreading of the 
front into the outer regions. However, the SP3 solution clearly has a smaller and sharper 
front than the other two lower-order approximations. Quantitatively, this becomes more 



clear when looking at a cut through the energy profile at x = 3.5, shown in Figure 4(b) 
The escaping front to the top is between y = 5 and y = 7 and it becomes clear that SP3 
is closest to the benchmark, with SSP3 having slight advantages over SPi. 



6 Conclusions 

Concerning the validity of the time-dependent SPn equations, assertions similar to the 
steady case can be made. Physically, the parabolic scaling and e small mean that we require 
the ratio of mean free path and characteristic length scale, as well as the characteristic 
length divided by the product of characteristic time and velocity to be small of order e. 

The numerical results here and in [7| indicate that the SPn approximations improve diffu- 
sion theory in the sense that not too far away from the diffusive limit a better approximation 
is obtained. 

In the lattice test case, a strong grid refinement around the obstacles was necessary. Sim- 
ilarly, to approximate the effect of source terms well, initially smaller time steps were 
needed. Using adaptivity in both time and space discretization is essential in many radia- 
tive transfer applications. This is especially true for non-homogeneous media with varying 
coefficients or with source terms. 
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